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Abstract 

We investigate the form and dynamics of shock-acoustic waves generated 
by earthquakes. We use the method for detecting and locating the sources of 
ionospheric impulsive disturbances, based on using data from a global network 
of receivers of the GPS navigation system and requiring no a priori information 
about the place and time of associated effects. The practical implementation 
of the method is illustrated by a case study of earthquake effects in Turkey 
(August 17, and November 12, 1999), in Southern Sumatera (June 4, 2000), 
and off the coast of Central America (January 13, 2001). It was found that 
in all instances the time period of the ionospheric response is 180-390 s, and 
the amplitude exceeds by a factor of two as a minimum the standard deviation 
of background fluctuations in total electron content in this range of periods 
under quiet and moderate geomagnetic conditions. The elevation of the wave 
vector varies through a range of 20-44°, and the phase velocity (1100-1300 
m/s) approaches the sound velocity at the heights of the ionospheric F-region 
maximum. The calculated (by neglecting refraction corrections) location of the 
source roughly corresponds to the earthquake epicenter. Our data are consistent 
with the present views that shock-acoustic waves are caused by a piston-like 
movement of the Earth surface in the zone of an earthquake epicenter. 
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1 Introduction 



A plethora of publications have been devoted to the study of the ionospheric response 
to disturbances arising from impulsive forcing on the Earth's atmosphere. It was 
found that in many cases a high proportion of the energy of the initial atmospheric 
disturbance is concentrated in the acoustic shock wave. Large earthquakes also 
provide a natural source of impulsive forcing. 

These investigations have also important practical implications since they furnish 
a means of substantiating reliable signal indications of technogenic effects (among 
them, unauthorized), which is necessary for the construction of an effective global ra- 
diophysical system for detection and localization of these effects. Essentially, existing 
global systems of such a purpose use different processing techniques for infrasound 
and seismic signals. However, in connection with the expansion of the geography and 
of the types of technogenic impact on the environment, very challenging problems 
heretofore have been those to improve the sensitivity of detection and the reliability 
of measured parameters of the sources of impacts, based also on using independent 
measurements of the entire spectrum of signals generated during such effects. 

To solve the above problems requires reliable information about fundamental 
parameters of the ionospheric response to the shock wave, such as the amplitude and 
the form, the period, the phase and group velocity of the wavetrain, as well as angular 
characteristics of the wave vector. Note that for naming the ionospheric response 
of the shock wave, the literature uses terminology incorporating a different physical 
interpretation, among them the term 'shock-acoustic wave' (SAW) (Nagorsky, 1998). 
For convenience, in this paper we shall use this term despite the fact that it does 
not reflect essentially the physical nature of the phenomenon. 

There is a wide scatter of published data on the main parameters of SAWs 
generated during industrial explosions and earthquakes. The oscillation period of 
the ionospheric response of SAWs varied from 30 to 300 sec, and the propagation 
velocity fluctuated from 700 to 1200 m/s (Nagorsky, 1998; Fitzgerald, 1997; Calais 
et al., 1998; Afraimovich et al., 1984; Blanc and Jacobson, 1989). 
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The lack of comprehensive, reliable data on SAW parameters is due primarily 
to the limitations of existing experimental methods and detection facilities. The 
main body of data was obtained by measuring the frequency Doppler shift at ver- 
tical and oblique-incidence ionospheric soundings in the HF range (Nagorsky, 1998; 
Afraimovich et al., 1984; Jacobson and Carlos, 1994). In some instances the sensi- 
tivity of this method is sufficient to detect the SAW reliably; however, difficulties 
emerge for localizing the region where the detected signal is generated. These prob- 
lems are caused by the multiple-hop character of HF signal propagation. This gives 
no way of deriving reliable information about the phase and group velocities of the 
SAW propagation, estimating angular characteristics of the wave vector, and, further 
still, of localizing the SAW source. 

Using the method of transionospheric sounding with VHF radio signals from 
geostationary satellites, a number of experimental data on SAW parameters were 
obtained in measurements of the Faraday rotation of the plane of signal polarization 
which is proportional to a total electron content (TEC) along the line connecting 
the satellite-borne transmitter with the receiver (Li et al., 1994). 

A common drawback of the above-mentioned methods when determining the 
SAW phase velocity is the necessity of knowing the time of events since this velocity 
is inferred from the SAW delay with respect to the time of events assuming that the 
velocity is constant along the propagation path, which is quite contrary to fact. 

For determining the above-mentioned rather complete set of SAW parameters, 
it is necessary to have appropriate spatial and temporal resolution which cannot 
be ensured by existing, very sparse, networks of ionosondes, oblique-incidence radio 
sounding paths, and incoherent scatter radars. 

The advent of the Global Positioning System (GPS) plus the subsequent creation 
of extensive networks of GPS stations (at least 757 sites as of November 2000), with 
their data being now available via the Internet, opened up a new era in remote 
ionospheric sensing. 
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Currently some authors have embarked on an intense development of methods 
for detecting the ionospheric response of strong earthquakes (Calais and Minster, 
1995), rocket launchings (Calais and Minster, 1996), and industrial surface explo- 
sions (Fitzgerald, 1997; Calais et al., 1998). In the cited references the SAW phase 
velocity was determined by the 'crossing' method by estimating the time delay of 
SAW arrival at subionospheric points corresponding to different GPS satellites ob- 
served at a given time. However, the accuracy of such a method is rather low 
because the altitude at which the subionospheric points are specified, is determined 
in a crude way. 

The goal of this paper is to describe a method for determining parameters of the 
SAW generated by earthquakes (including the phase velocity, angular characteristics 
of the SAW wave vector, the direction towards the source, and the source location) 
using GPS-arrays whose elements can be chosen out of a large set of GPS stations 
from the global GPS network. 

Section 2 presents a description of the experimental geometry, and general infor- 
mation about the earthquakes under consideration. The proposed method is briefly 
outlined in Section 3. Results of measurements of SAW parameters from different 
GPS arrays during earthquakes are presented in Section 4. Section 5 is devoted to 
the discussion of experimental results, including analytical simulation results. 



2 The geometry and general characterization of exper- 
iments 

Detection results on two earthquakes in Turkey (August 17, and November 12, 
1999), in Southern Sumatera (June 4, 2000), and off the coast of Central America 
(January 13, 2001) are presented below. The information about the earthquakes, 
was acquired via the Internet ' http:/ /earthquake .usgs.gov| '. General information 



about these earthquakes is presented in Table 1 (including the time of the main 
shock to in the universal time UT, the position of the earthquake epicenter, depth, 
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the magnitude, as well as the level of geomagnetic disturbance from the data on 
Dst- variations). It was found that the deviation of Dst for the selected days was 
quite moderate, which enabled the SAWs to be identified. 

Fig. 1 illustrates the experimental geometry during the earthquakes in Turkey - 
a) , off the coast of Central America - b) , and in Southern Sumatera - c) . 

In spite of the small number of GPS stations in the earthquake area, we were able 
to use a sufficient number of them for the implementation of the proposed method. 

Table 2 presents the geographic coordinates of the GPS stations used as GPS 
array elements. 

3 Methods of determining shock-acoustic wave charac- 
teristics using GPS-arrays 

The standard GPS technology provides a means for wave disturbances detecion based 
on phase measurements of TEC at each of spaced two-frequency GPS receivers. A 
method of reconstructing TEC variations was detailed and validated in a series of 
publications (Calais and Minster, 1995, 1996; Fitzgerald, 1997). We reproduce here 
only the final formula for the total electron content (I) 

1 = 40^08 W^fl [{LlXl ~ L " X2) + C ° nSt + UL] (1) 
where LiAi and L2X2 are additional paths of the radio signal caused by the phase 
delay in the ionosphere, (m); L\ and Li represent the number of phase rotations at 
the frequencies f\ and fi, \\ and A2 stand for the corresponding wavelengths, (m); 
const is the unknown initial phase ambiguity, (m) ; and nL are errors in determining 
the phase path, (m). 

Phase measurements in the GPS can be made with a high degree of accuracy 
corresponding to the error of TEC determination of at least 10 14 m -2 when averaged 
on a 30-second interval, with some uncertainty of the initial value of TEC, however. 
This makes possible detection of ionization irregularities and wave processes in the 
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ionosphere over a wide range of amplitudes (up to 10 -4 of the diurnal TEC variation) 
and periods (from 24 hours to 5 min). The unit of TEC, TECU, is equal to 10 16 m -2 , 
it is commonly accepted in the literature, and will be used in the following. 

A convenient way of determining the ionospheric response delay of the shock wave 
involves the frequency Doppler shift F from TEC series obtained by formula (|l|). 
Such an approach is also useful in comparing TEC response characteristics from 
the GPS data with those obtained by analyzing VHF signals from geostationary 
satellites, as well as in detecting the shock wave in the HF range. To an approxima- 
tion sufficient for the purpose of our investigation, a corresponding relationship was 
obtained by K. Davies (1969) 

F = 13.5 • l(T 8 / t 7/ (2) 

where 1[ stands for the time derivative of TEC. Relevant results derived from an- 
alyzing the .F(i)-variations calculated for the 'reduced' frequency of 136 MHz are 
discussed in Section |I[ 

The correspondence of space-time phase characteristics, obtained through tran- 
sionospheric soundings, with local characteristics of disturbances in the ionosphere 
was considered in detail in a wide variety of publications (Afraimovich et al., 1992; 
Mercier and Jacobson, 1997) and is not analyzed at length in this study. The most 
important conclusion of the cited references is the fact that, as for the extensively 
exploited model of a 'plane phase screen' disturbances AI(x, y, t) of TEC faithfully 
copy the horizontal part of the corresponding disturbance of local electron concen- 
tration AN(x,y, z,t), independently of the angular position of the source, and can 
be used in experiments for measuring the TEC disturbances. 

However, the TEC response amplitude experiences a strong azimuthal depen- 
dence caused by the integral character of a transionospheric sounding. As a first 
approximation, the transionospheric sounding method is responsive only to Travel- 
ing Ionospheric Disturbances (TIDs) with the wave vector Kt perpendicular to the 
direction r wich is along the Line-of-Sight (LOS) from the receiver to the satellite. 
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A corresponding condition for elevation 9 and azimuth a of an arbitrary wave vector 
Kt normal to the direction r, has the form 

9 = arctan(— cos(a s — a)/ tan 9 S ) (3) 

where a s is the azimuthal angle measured east to north, and 9 S is the angle of 
elevation of the satellite at the receiver. 

We used formula (|3|) to determine the elevation 9 of Kt from the known mean 
value of azimuth a by Afraimovich et al. (1998) - see Sections |3.2| and ||. 

3.1 Detection and determination of the horizontal phase velocity 
and the direction of the SAW phase front along the ground by 
GPS- arrays 

In the simplest form, space-time variations of the TEC A/(i, x, y) in the ionosphere, 
at each given time t can be represented in terms of the phase interference pattern 
that moves without a change in its shape (the solitary, plane travelling wave) 

AI(t, x,y) = 5 sm(Qt - K x x - K y y + ip ) (4) 

where 5, K x , K y , Q, are the amplitude, the x- and y-projections of the wave vector 
K, and the angular frequency of the disturbance, respectively; T = 2tt/VL, A = 
2n/\K\ is its period and wavelength; and ipo is the initial phase of the disturbance. 
The vector K is a horizontal projection of the full vector Kt. 

At this point, it is assumed that in the case of small spatial and temporal in- 
crements (the distances between GPS-array sites are less than the typical spatial 
scale A of TEC variation, and the time interval between counts is less than the 
corresponding time scale), the influence of second derivatives can be neglected. The 
following choices of GPS-arrays all meet these requirements. 

We now summarize briefly the sequence of data processing procedures. Out of 
a large number of GPS stations, three sites (A, B, C) are selected, within distances 
not exceeding about one-half of the expected wavelength A of the perturbation. Site 
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B is taken to be the center of a topocentric reference frame whose axis x is directed 
east, and the axis y is directed north. The receivers in this frame of reference have 
the coordinates (xa, Va), (0,0), (xc, yc)- Such a configuration of the GPS receivers 
represents the GPS-array with a minimum number of the required elements. In 
regions with a dense network of GPS sites, we can obtain a large variety of GPS- 
arrays of a different configuration enabling the acquired data to be checked for 
reliability; in this paper we have exploited this possibility. 

The input data include series of slant TEC values /^(t), I Bit)-, lit), as well as 
corresponding series of elevation values s (t) and the azimuth a s {t) of the LOS. For 
determining SAW characteristics, continuous series of measurements of /^(i), /b(£), 
Ic(t) are selected with a length of at least a one-hour interval which includes the 
time of a earthquake. 

To eliminate spatio-temporal variations of the regular ionosphere, as well as 
trends introduced by orbital motion of the satellite, a procedure is used to remove 
the trend involving a preliminary smoothing of the initial series with the selected 
time window. This procedure is better suited to the detection of a single pulse 
signal (iV-wave) than the frequently used band-pass filter (Li et al., 1994; Calais 
and Minster, 1995, 1996; Fitzgerald, 1997; Calais et al., 1998). A limitation of 
the band-pass filter is the oscillatory character of the response which prevets from 
reconstructing the form of the N-w&ve. 

Elevation s (t) and azimuth a s (t) values of the LOS are used to determine the 
location of the subionospheric point, as well as to calculate the elevation 9 of the 
wave vector Kt of the disturbance from the known azimuth a (see formula (||)). 

The most reliable results from the determination of SAW parameters correspond 
to high values of elevations 6 s (t) of the LOS because sphericity effects become rea- 
sonably small. In addition, there is no need to convert the slant TEC AI(i) to a 
'vertical' value. In this paper, all results were obtained for elevations 9 s (t) larger 
than 30°. 
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Since the distance between GPS-array elements (from several tens of kilometers 
to a few hundred kilometers) is much smaller than that to the GPS satellite (over 
20000 km), the array geometry at the height of the ionosphere is identical to that 
on the ground. 

Fig. 2a shows typical time dependencies of an slant TEC I(t) at the GPS-array 
BSHM station near the area of the earthquake of August 17, 1999. (heavy curve), 
one day before and after the earthquake (thin lines). For the same days; panel 
b shows TEC variations AI(t) after removal of a linear trend and smoothing by 
averaging over a sliding window of 5-min. Variations in frequency Doppler shift 
F(t), 'reduced' to the sounding signal frequency of 136 MHz, for three sites of the 
array (KATZ BSHM GILB) on August 17, 1999, are presented in panel c. 

Fig. 2 shows that fast iV-shaped oscillations, with a typical period of about 
390 s, are distinguished among slow TEC variations. The oscillation amplitude (up 
to 0.12 TECU) is far in excess of the background TEC fluctuation intensity as seen 
on the days before and after the earthquake. Variations in frequency Doppler shift 
F(t) for spatially separated sites (KATZ BSHM GILB) are well correlated but are 
shifted relative to each other by an amount well below the period, which permits the 
SAW propagation velocity to be unambiguously determined. The 30 s sampling rate 
of the GPS data is not quite sufficient for determining small shifts of such signals 
with an adequate accuracy for different sites of the array. Therefore, we used a 
parabolic approximation of the .F(i)-oscillations in the neighborhood of minimum 
F(t), which is quite acceptable when the signal/noise ratio is high. 

Taking into account the good signal/noise ratio (better than 1), and knowing the 
coordinates of the array sites A, B and C, we determine the horizontal projection of 
the phase velocity Vh from time shifts t p of a maximum deviation of the frequency 
Doppler shift F(t). Preliminarily measured shifts are subjected to a linear transfor- 
mation with the purpose of calculating shifts for sites spaced relative to the central 
site northward N and eastward E. This is followed by a calculation of the E- and 
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iV-components of V x and V y , as well as the direction a in the range of angles 0-360° 
and the modulus V/ 4 of the horizontal component of the SAW phase velocity 

a = arctan(V^/V^) 

(5) 

V h =\V x V y \(V* + V*)- 1 / 2 
where V y , V x are the velocities with which the phase front crosses the axes x and y. 
The orientation a of the wave vector K, which is coincident with the propagation 
azimuth of the SAW phase front, is calculated unambiguously in the range 0-360° 
subject to the condition that arctan{V y /V x ) is calculated having regard to the sign 
of the numerator and denominator. 

The above method for determining the SAW phase velocity neglects the correc- 
tion for orbital motion of the satellite because the estimates of Vh obtained below 
exceed an order of magnitude as a minimum the velocity of the subionospheric point 
at the height of the ionosphere for elevations 6 S > 30° (Afraimovich et al., 1998). 

Obviously the method presented above can be used if the distance between GPS 
stations is much shorter than the TEC disturbance wavelength A and the distance 
from the earthquake epicenter to the array. This corresponds to the detection con- 
dition in the far-field zone. 

From the delay At = t p — to and the known path length between the earthquake 
focus and the subionospheric point we calculated also the SAW mean velocity V a 
in order to compare our estimates of the SAW phase velocity with the usually used 
method of measuring this quantity. 

3.2 Determination of the wave vector elevation 9 and the velocity 
modulus V t of the shock wave 

Afraimovich et al. (1992) showed that for the Gaussian ionization distribution the 
TEC disturbance amplitude (M) is determined by the aspect angle 7 between the 
vectors K t and r, as well as by the ratio of the wavelength of the disturbance A to 
the half-thickness of the ionization maximum 
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M oc exp 



( 



Tr 2 h d cos 2 7 



A 2 cos 2 9 



) 



(6) 



In the case under consideration (see below), for a phase velocity on the order 
of 1 km/s and for a period of about 200 s, the wavelength A is comparable with 
the half-thickness of the ionization maximum h d . When the elevations 9 S are 30°, 
45°, 60° , the 'beam- width' M( 7 ) at 0.5 level is, respectively, 25°, 22° and 15°. If 
hd is twice as large as the wavelength, then the beam tapers to 14°, 10° and 8°, 
respectively. 

The beam-width is sufficiently small that the aspect condition (|3|) restricts the 
number of beam trajectories to the satellite, for which it is possible to detect reliably 
the SAW response in the presence of noise (near the angles 7 = 90° ) . On the other 
hand, formula (||) can be used to determine the elevation 9 of the wave vector Kt 
of the shock wave at the known value of the azimuth a (Afraimovich et al., 1998). 
Hence the phase velocity modulus Vt can be defined as 



The above values of the width M(-f) determine the error of calculation of the 
elevations 6 (of order 20° to the above conditions). 

3.3 Determining the position of the SAW source without regard 
for refraction corrections 

The ionospheric region that is responsible for the main contribution to TEC varia- 
tions lies in the neighborhood of the maximum of the ionospheric F-region, which 
does determine the height hi of the penetration point. When selecting hi, it should 
be taken into consideration that the decrease in electron density with height above 
the main maximum of the i<2-layer proceeds much slower than below this maxi- 
mum. Since the density distribution with height is essentially a 'weight function' of 
the TEC response to a wave disturbance (Afraimovich et al., 1992), it is appropriate 
to use, as hi, the value exceeding the true height of the layer h m F2 maximum by 



V t = V h - cos(9) 



(7) 
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about 100 km. h m F2 varies between 250 and 350 km depending on the time of day 
and on some geophysical factors which, when necessary, can be taken into account 
if additional experimental data and current ionospheric models are available. In all 
calculations that follow, hi = 400 km is used. 

To a first approximation, it can be assumed that the imaginary detector, which 
records the ionospheric SAW response in TEC variations is located at this altitude. 
The horizontal extent of the detection region, which can be inferred from the prop- 
agation velocity of the subionospheric point as a consequence of the orbital motion 
of the GPS satellite (on the order of 70 -150 m/s; see Pi at al., 1997), and from the 
SAW period (on the order of 200 s — see Section does not exceed 20-40 km, 
which is far smaller than it 'vertical size' (on the order of the half-thickness of the 
ionization maximum h^). 

From the GPS data we can determine the coordinates X s and Y s of the subiono- 
spheric point in the horizontal plane X0Y of a topocentric frame of reference centered 
on the point B(0,0) at the time of a maximum TEC deviation caused by the arrival 
of the SAW at this point. Since we know the angular coordinates 9 and a of the wave 
vector Kt , it is possible to determine the location of the point at which this vector 
intersects the horizontal plane X'OY' at the height h w of the assumed source. As- 
suming a rectilinear propagation of the SAW from the source to the subionospheric 
point and neglecting the Earth sphericity, the coordinates X w and Y w of the source 
in a topocentric frame of reference can be defined as 

X w = X p — (hi — h w ) : — - — (8) 

Y w = Y p -( hi -h w )^^ (9) 

smu 

The coordinates X w and Y w , thus obtained, can readily be recalculated to the 
values of the latitude and longitude (<f> w and X w ) of the source. 

For SAW generated during earthquakes, industrial explosions and underground 
tests of nuclear devices, h w is taken to be equal to (the source lying at the ground 
level) . 
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4 Results of measurements 

Hence, using the transformations described in Section ||, we obtain the parameters 
set determined from TEC variations and characterizing the SAW (see Table 3). 

Let us consider the results derived from analyzing the ionospheric effect of SAW 
during earthquake August 17, 1999 obtained at the array (KATZ, BSHM, GILB) 
for PRN 6 (at the left of Fig. 2, and line 2 in Table 3). 

In this case the delay of the SAW response with respect to the time of the 
earthquake is 20 min (DAY 229). The SAW has the form of an iV-wave with a 
period T of about 390 s and an amplitude Aj = 0.12 TECU, which is an order of 
magnitude larger than TEC fluctuations for background days (DAY 228, DAY 230). 
It should be noted, however, that this time interval was characterized by a very low 
level of geomagnetic activity (-14 nT). 

The amplitude of a maximum frequency Doppler shift Ap at the 'reduced' fre- 
quency of 136 MHz was found to be 0.04 Hz. In view of the fact that the shift F is 
inversely proportional to the sounding frequency squared (Davies, 1969), this corre- 
sponds to a Doppler shift at the working frequency of 13.6 MHz and the equivalent 
oblique-incidence sounding path of about Ap = 4 Hz. 

Solid curves in Fig. la represent the trajectories of the subionospheric points for 
each GPS satellite at the height /ij=400 km during the time interval 0.0-0.8 UT 
for August 17, 1999, and 17.0-17.4 UT for November 12, 1999. Dark diamonds 
along the trajectories correspond to the coordinates of the subionospheric points at 
time t p of a maximum deviation of the TEC (Fig. 2b and Fig. 2e). Crosses show 
the positions of the earthquake epicenters. Asterisks mark the source location at 
km altitude inferred from the data from the GPS arrays. Numbers at the asterisks 
correspond to the respective day numbers. Straight dashed lines that connect the 
expected source and the subionospheric point represent the horizontal projection of 
the wave vector Kf. 

The azimuth a and elevation 9 of the wave vector Kt whose horizontal projec- 
tion is shown in Fig. la by a dashed line and is marked by K%, are 154° and 26°, 
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respectively. The horizontal component and the modulus of the phase velocity were 
found to be Vh = 1307 m/s and Vt = 1174 m/s. The source coordinates at km 
altitude were determined as <p w = 39.1° and X w = 25.9°. The alculated (by neglect- 
ing refraction corrections) location of the source roughly was close the earthquake 
epicenter. 

The 'mean' velocity of about V a = 870 m/s, determined in a usual manner from 
the response delay with respect to the start, was smaller then the phase velocity Vt. 
Conceivably this is associated with an added delay of the response as a consequence 
of the refraction distortions of the SAW path along the LOS which were neglected 
in this study. 

Similar results for the array (KABR, TELA, GILB) and PRN 30 were also ob- 
tained for the earthquake of November 12, 1999. They correspond to the projection 
of the vector K<i in Fig. la, the time dependencies in Fig. 2 at the right, and line 
7 in Table 3. The only point deserving mention here is that the SAW amplitude 
was somewhat smaller than that for the earthquake of August 17, 1999. With an 
increased level of geomagnetic activity ( -44 nT), this led to a smaller (compared 
with August 17, 1999) signal/noise ratio; yet this did not preclude reliable estimates 
of the SAW parameters. 

A comparison of the data for both earthquakes showed a reasonably close agree- 
ment of SAW parameters irrespective of the level of geomagnetic disturbance, the 
season, and the local time. 

To convince ourselves that the determination of the main parameters of the 
SAW form and dynamics is reliable for the earthquakes analyzed here, in the area 
of the earthquake we selected different combinations of three sites out of the sets of 
GPS stations available to us, and these data were processed with the same processing 
parameters. Relevant results (including the average results for the sets £), presented 
in Table 3 and in Fig. la (SAW source position), show that the values of SAW 
parameters are similar, which indicates a good stability of the data, irrespective of 
the GPS-array configuration. 
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The aspect condition (||), corresponding to a maximum amplitude of the TEC 
response to the transmission of SAWs, was satisfied quite well for this geometry, 
simultaneously for all stations. This is confirmed by a high degree of correlation 
of the SAW responses at the array elements (Fig. 2c and Fig. 2f), which made it 
possible to obtain different sets of triangles out of the six GPS stations available to 
us. 

The relative position of the GPS stations was highly convenient for determining 
the SAW parameters during the earthquakes in Turkey, and met the implementation 



conditions for the method described in Section |3.1| . Thus, the distance between 
stations (100-150 km at most) did not exceed the SAW wavelength of about 200-300 
km, and was far less than the distance from the epicenter to the array (1000 km). 

Let us consider the results derived from analyzing the ionospheric effect of SAW 
during earthquake January 13, 2001, obtained at the array (TEGU, MANA, ESTI) 
for PRN13 (at the left of Fig. 3 and line 11 in Table 3). 

In this case the delay of the SAW response with respect to the time of the 
earthquake is 15 min (DAY 013). The SAW has the form of oscillations with a 
period T of about 270 s and an amplitude Aj = 0.2 TECU, which is an order of 
magnitude larger than TEC fluctuations for background days (DAY 012, DAY 014). 
It should be noted, that this time interval was characterized by a very low level of 
geomagnetic activity (4 nT). 

Similar results for the array (SAMP, NTUS, BAKO) and PRN 03 were also 
obtained for the earthquake of June 4, 2000. (at the right of Fig. 3 and line 10 in 
Table 3). 

Note that in this case the response amplitude exceeded twice, as a minimum, 
that for the other events under consideration. It is not improbable that this is due 
to the maximum magnitude of the earthquake in Southern Sumatera (see Table 1). 

Unfortunately, because of the inadequately well-developed network of stations, 
for the earthquakes in Southern Sumatera (June 4, 2000) and off the coast of Central 
America (January 13, 2001) for each of these events it was impossible to select arrays 
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meeting the applicability conditions of the method for determining the SAW wave 



vector parameters described in Section |3.1 . 

It is evident from the geometry in Fig. lb and Fig. lc that the earthquake epi- 
centers lay inside the GPS arrays, which does not meet the far-field zone condition. 
It is possible that this is also responsible for the form of the response (presence of 
strong oscillations) which differs from the N-form of the response for the earthquakes 
in Turkey. On this basis, for these earthquakes we can point out only the very fact 
of reliable detection of the response, and determine its amplitude Aj, typical period 
T, delay At = t p — to, and velocity V a (see Table 3). The values of these quantities 
were close to the data obtained for the earthquakes in Turkey. 



5 Discussion 

The data of Table 3 are quite sufficient to estimate the position of the disturbances 
under discussion in the diagnostic diagram of the atmospheric waves. Specifically, 
the values of the characteristic periods of bipolar signals, presented in Fig. 2b and 
Fig. 2e, are T\ ~ 300 s, and T2 ~ 200 s. The wave vectors are at the angles 0% ~ 70° 
and /?2 ~ 50° with vertical. At the height z ~ 400 km, in turn, the periods, 
corresponding to local frequencies of the acoustic cut-off oj a and the Brunt- Vaisala 
u>b, are: T a = 2i:/oj a « 950 s, T& = 2tt /uj^ rj 1050 s. 

In an isothermal atmosphere, only harmonics with periods larger than T^j sin (5 
can be assigned to the branch of internal gravity waves. This value exceeds sig- 
nificantly the values of T\ 2- Furthermore, is considerably smaller than T a . 
Hence we can contend that under conditions of the real (nonisothermal) atmosphere 
the disturbances under discussion almost entirely pertain to the branch of acoustic 
waves. 

Since in (short period) acoustic waves the variation AN in neutral density N 
satisfies the relation AN/N ~ u/C (u is the gas velocity in the wave, and C is the 
local velocity of sound), it is possible to make a lower estimate of the intensity u/C 
of the waves at the expected altitudes z e ff ~ 350-400 km: AN/N as AN e /N e ~ 
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AI/I « 0.02 — 0.04. In this estimation, the influence of the magnetic field is omitted, 
and the index e refers to electrons. The actual intensity of the acoustic wave must 
be higher than the estimated value presented above. 

In the case of an earthquake, the movements of the terrestrial surface are plau- 
sible sources of acoustic waves. A generally known source of the first type is the 
Rayleigh surface wave propagating from the epicentral zone. The ground motion in 
the epicentral area is the source of the second type. Let us discuss these possibilities. 

5.1 The Rayleigh wave 

The phase propagation velocity of the Rayleigh wave Vr « 3.3 km/s. Since Vr S> Co, 
where Co ~ 0.34 km/s is the sound velocity at the ground, only acoustic waves can be 
emitted (Golitsyn and Klyatskin, 1967). At a sufficient distance from the epicenter 
where the curvature of the Rayleigh wave front can be neglected, and with the proviso 
that uj 3> u a , acoustic waves are emitted at the angle (3r arcsin(Co/VR) « 6° to 
the vertical. 

Rayleigh waves propagate generally in the form of a train consisting of several 
oscillations whose typical period only rarely exceeds several tens of seconds. Acoustic 
waves are emitted upward in the form of the same train. Because of a strong 
absorption of the periodic wave, the only thing that is left over in the case of the 
acoustic train at heights z > 350 km is the leading phase of compression. It seems 
likely that only in the case of strong earthquakes (Alaskian earthquake of 1964), 
even at large distances from the epicenter, the disturbance (the leading portion of 
the acoustic train) that remains from the acoustic train, can have at these altitudes 
a duration of about 100 s, and quite an appreciable (u/C > 0.1) intensity (Orlov 
and Uralov, 1987) the nonlinear acoustics approximation. 

Unfortunately, the parameters of Rayleigh waves in the neighborhood of the 
subionospheric points appearing in Table 3 are unknown to us. However, a most pro- 
nounced bipolar character of the main signal (Fig. 2b, Fig. 2e, Fig. 3b and Fig. 3e), 
its intensity, and its long duration cast some doubt upon the fact that it is the 
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Rayleigh wave which is responsible for its origin. Such a conclusion is consistent 
with the observed propagation direction of the bipolar pulse, which makes a large 
angle with the vertical: (3 1>2 ~ 70-50° > R « 15-18°. 

Here it is taken into consideration that the acoustic ray originating from a point 
on the ground at the angle of (3r « 6° now makes at the heights z = z e ff an 
angle (3" R > (3r because of the refraction effect in a standard atmospheric model: 
Co/sin (3r = C /sin (3" R = const = Vr, C{z = z e ff) ~ 0.9-1 km/s. 

The presence of strong winds at ionospheric heights can alter the value of (3r. 
However, irrespective of the atmospheric model, the phase velocity Vh (Table 3) 
of the horizontal trace of the acoustic disturbance generated by the Rayleigh wave 
must coincide with Vr, and this is also not observed: Vh <C Vr. 

5.2 The epicentral emitter 

The detection of ionospheric disturbances, which are presumably generated by a 
vertical displacement of the terrestrial surface directly in the epicentral zone of an 
earthquake, using the GPS probing method, is reported in (Calais and Minster, 
1995). 

Results of the present study lend support to the above conjecture. However, the 
specific formation mechanism for the disturbance itself is still unclear. An approach 
to solving this problem is contained in earlier work, and involves substituting the 
epicentral emitter for a surface velocity point source or an explosion. In particular, 
the substitution of the earthquake zone for a point source turns out to be fruit- 
ful when describing long-period internal gravity waves at a very long (thousands 
of kilometers) distance from the epicenter (Row, 1967). The visual resemblance of 
ionospheric disturbances at short (hundreds of kilometers) distances from the earth- 
quake epicenter to disturbances from surface explosions is discussed in (Calais et al., 
1998). 

It should be noted that ionospheric disturbances generated by industrial surface 
and underground nuclear explosions are also visually similar. However, the genera- 
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tion mechanisms for disturbances are fundamentally different in this case (Rudenko 
and Uralov, 1995). The radiation source in underground nuclear tests is, as in the 
case of earthquakes, the terrestrial surface disturbed by the explosion. The intensity 
and spectral composition of the generated acoustic signal reveal a strong (unlike the 
surface explosion) dependence on the zenith angle, and are wholly determined by 
the form, size and characteristics of the movement of the terrestrial surface in the 
epicentral zone of the underground explosion. 

In this Section, we shall propose a model which, we hope, will help to understand 
the generation mechanism for acoustic disturbances — the subject of this paper. 
Because of the complexity of the problem formulated, and for lack of sufficient data 
on characteristics of the movement of the terrestrial surface in epicentral zones of 
earthquakes, the idealized model under discussion has an illustrative character. The 
computational scheme proposed below represents a simplified variant of the scheme 
used in Rudenko and Uralov, (1995) to calculate ionospheric disturbances generated 
by an underground confined nuclear explosion. 

5.3 The problem of radiation of the acoustic signal 

For the sake of simplicity, we consider a problem having an axial symmetry about 
the vertical axis z passing through the earthquake epicenter r = z = 0. The epi- 
central emitter is a set of plane annular velocity sources with the specified law of 
motion along the vertical U(r,t). Since our interest is with the estimation of the 
characteristics of an acoustic disturbance at a sufficient distance from the emitter, 
we take advantage of the far-field approximation of a linear problem of radiation. In 
the approximation of linear acoustics to 0J a and with no absorption present, the 
gas velocity profile in the wave can be estimated by the expression: 



where r = t — J^dl/C, y = rs ™P . Here (5 is the zenith angle of departure of the 
acoustic ray from the point r = z = 0. I is the group path length (the distance 




(10) 
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along the ray). a(t ,r) = ^jf'^ is the vertical acceleration of the terrestrial surface. 



A = A(z) = J%2o j s aC oustic factor, po, p stand for the air density at the 
ground and at the height z, respectively, and L is a typical radius of the epicentral 
emitter. 

In an isothermal atmosphere where the ray trajectories are straight lines, the 
quantity R = (V r 2 + z 2 = I) is the radius of a divergent (in the far-field approxima- 
tion) spherical wave. In this case the cross-section of the selected ray tube S oc R 2 . In 
the real atmosphere the value of S is determined from geometrical optics equations, 
and using the expression ( |l0| ) requires a further complication of the computational 
scheme. Nevertheless, to make estimates we shall use only the relation (|To|) , and, in 
doing so, the value of R will be corrected. 

The case = (and also I = R = z) is a special one: 



A f L 

u(t, = 0) = — — / a(r, r)rdr (11) 



RCqjq 

Let us discuss the situation where all ring-type emitters 'operate' synchronously: 
U(r,t) = U{t). In this case the epicentral zone is emitting as a round piston 2L in 
diameter. Assume also that with a shock of the earthquake, the vertical displacement 
£, the velocity U = d^/dt and the acceleration a = dU/dt are time dependent, as 
shown in Fig. 4a. For a rectangular velocity impulse, i.e. in the limit At 5t — > 
(in this case a5t — ► ±[/o<5(f'), where 5(t) is a delta-function), from (10) we can 
obtain 



U C A 



M(T ' /3) = vrW^ ^i " r2||T| ^ " V f£ - ( T - At) 2 || T -At|< y J (12) 
where = Ls ^f \ Uo, At are the amplitude and duration of the rectangular velocity 
impulse. In this case the vertical displacement of the terrestrial surface after the 
earthquake shock is £o = UoAt. In the strict sense, the expression (12) holds true 
if the condition 3> 5t is satisfied, i.e. with zenith angles 3> 0* ~ arcsin^j^. 
Here 8t is the operation duration of the terrestrial shock at the beginning of the 
movement and at the stop of the piston. When w 0* <C 1 , the expression (|i~2|) and 
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the expression (11), in view of 5t ^ 0, yield approximately identical signals. What 



actually happens is that the validity of the expression ( |12[) will break down still 
earlier, and we are justified in using it only for zenith angles j3 > f3** = arctan(^) > 
13*. 

The curves in Fig. 4b give an idea of the relative amplitude and form of acoustic 
signals u(t, 0) in the far-field zone of radiation of the piston 2L = 60 km in diameter. 
The duration of the rectangular velocity impulse of the piston was chosen arbitrarily, 
At = 10 s. The signals correspond to the expression (^) under the assumption of 
a homogeneous, A = 1, C = Co ~ 0.34 km/s, atmosphere. The zenith angles are 
(3 = 10°, 20°, 40°, 60°, 90°. The spherical surface R = const serves as a reference. 
The abscissa axis indicates the time r in seconds, and the axis of ordinates indicates 
the gas velocity u. In this case the amplitude of the most intense signal (/? = 10°) 
is taken to be unity. 

With increasing /3, the amplitude of the signals Umax 

decreases, and the duration 

T increases. The leading and trailing edges of the bipolar pulses in Fig. 2 and Fig. 3 
have an equal duration At. The positive part of the bipolar pulse corresponds to 
the compression phase of the acoustic wave, and the negative part refers to the 
rarefaction phase. The area of the compression phase (in coordinates u, t) equals 
the largest displacement \+ °f a un h volume of the atmosphere in the direction of 
propagation (along the ray) of the wave. A total area of a bipolar pulse is zero: 
X+ = — X-- For acoustic signals, described by the expression (|i~2j), the following 
useful relations hold true: 



2AL 



Umax = u ^^R S i^~^ T] ~ ri2 '' T = 2 yL + At ( 13 ) 



AL / 1 At 

X+ = ^2^R^ { y l - 11+ -r, arCSm ^ ^=2^ (14) 
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5.4 The problem of acoustic signal propagation 

The wave vectors Kt of the bipolar pulses make with the vertical the angles (3\ = 70° 
and = 50°. With the adopted values of z e ff = 350-400 km, the distances of 
the corresponding subionospheric points from the earthquake epicenters in Fig. la 
are approximately r\ = 800 km and T2 = 600 km. The rays, constructed in the 
approximation of linear geometrical acoustics (LGA) and having at the heights z = 
z e ff (C(z = z e ff) ~ 0.9-1 km/s) the propagation angles 0% and (3%, correspond to 
the zenith angle of departure ~ 19-21° and (3 ~ 15-17° at the level z = 0. 

The fact that these values of (3 satisfy the inequality (3 < (3* ~ 25-30°, is in 
reasonably good agreement with the familiar picture of rays from a ground-level 
point source — the rays with (3 > (3* are captured by the atmospheric waveguide 
z < z* = 120 km, and only the rays emitted upward inside the solid angle Q* ~ 
1 sterad can penetrate to the heights z > z* . In standard models of the atmosphere, 
however, to the values of the angles {3\p,{z = z e ff) there correspond the locations 
of subionospheric points lying several hundreds kilometers nearer to the epicenter 
compared with the experimental values of t\ ~ 800 km and T2 ~ 600 km. This 
inconsistency can be due to two reasons. 

One reason is that the value of z e ff ~ 350-400 km, which we are using, is too 
low. This is supported by the detection of velocities Vt ~ 1.2 km/s of traveling 
disturbances (Table 3) which exceed markedly the sound velocity C ~ 0.9-1 km/s 
at the heights of ~ 350-400 km. However, it seems likely that such a discrepancy 
may be disregarded, in view of the errors of the measurement technique used (the 
probability of an additional heating of the upper atmosphere prior to the earthquake 
cannot be ruled out, however). 

The increase of the actual value of z e ff can also be associated with a strong de- 
pendence (see (|i~3|), di~4] ) and Fig. 4b) of the power of the emitted signal on the zenith 
angle of departure of the ray from the earthquake epicenter. Verifying this factor 
requires a more detailed analysis, based on particular data on vertical movements of 
the terrestrial surface in the epicentral zone. Such data are unavailable to us. Below 
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is a discussion only of the second reason for the above-mentioned inconsistency as a 
highly probable one. 

The second reason may be associated with the violation of the validity conditions 
of the LGA approximation at a sufficient distance from the source of the acoustic 
disturbances under discussion. Indeed, the utilization of this approximation is justi- 
fied until the parameters of the medium and of the wave itself change substantially 
on the size of the first Fresnel zone dp ~ vAZ that determines the physical (trans- 
verse) size of the ray. Typical wavelengths of the bipolar pulses under discussion at 
the heights z = z e ff are large: A ~ 300-200 km. Distance I along the expected ray 
is of order 900-700 km. Then dp ~ 520-370 km, which exceeds substantially the 
scales of variation of atmospheric parameters. The value of dp is actually somewhat 
smaller, because the typical scale of a disturbance decreases as it approaches the 
source. A violation of the LGA approximation at the above-mentioned distances 
from the epicentral source occurs also for model signals [3 > 20°, shown in Fig. 4b. 

The increasing violation of the applicability conditions of the LGA approximation 
with an increase of I implies the transport of the wave energy not strictly along the 
calculated rays but along the lines with a smaller curvature. With a mere estimate of 
the dilution factor R in the expressions dH), (U) and (II), in mind, these lines will 
be assumed to be straight when z > z* and to originate from an imaginary source 
lying at the height z* ~ 120 km above the earthquake epicenter. Such a situation 
is also clearly manifested in the LGA approximation. At the heights z > z*, at 
t ~ 600 s, for example, the surface of the wave front from the ground-level impulsive 
source resembles the surface of a hemisphere centered on the point r ~ 0, z ~ z* . 

Ultimately the energy that arrives from below inside the solid angle O* ~ 1 is 
scattered into a solid angle O, ~ 2ir. From the condition of conservation of wave 
energy, it is possible to find Ri } 2 = \f^\/ (Az) 2 + r\ 2 , when Az = z e ff — z* . 
When using the expressions (|T3|), flT4| ) in the subsequent discussion, we will take the 
quantity R± rather than R, and the value of the zenith angle of departure will be 
taken to be j3 = 20°. These parameters approximately correspond to the generation 
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and propagation conditions of the signal from the first earthquake. The typical size 
of the epicentral zone of this earthquake was about 2L = 60 km (according to the 
USGS data: www.neic.cr.usgs.gov). This same value of L was used in calculating 
the signals shown in Fig. 4b. 

As is evident even from Fig. 4b (thick line), the duration of the signal, having 
a zenith angle of departure (3 = 20°, is about 70 s. When the signal propagates in 
the approximation of linear acoustics and with no absorption, its form and duration 
remain unchanging, and only its amplitude changes. 

In actual conditions the combined effect of the nonlinear attenuation and linear 
absorption factors leads to a stretching of the bipolar pulse, and to a change of its 
form (we do not discuss the dispersion factor). In this case the effect of the nonlin- 



early factor occurs in such a manner that the integral value of x+, (14), calculated 
as an approximation of linear acoustics, remains as such in the approximation of 
nonlinear acoustic as well. Moreover, taking into account the finite width AT s h of 
the shock front does not change this situation until AT s /j < T/4. The linear absorp- 
tion factor, in turn, somewhat reduces the true value of x+ because of the mutual 
diffusion of the compression and rarefaction phases. Nevertheless, for a hypothetical 
estimation of the earthquake parameters, we shall use the assumption about the 
conservation of the value of x+- 

As is intimated by Table 3, the mean value of the TEC disturbance amplitude 
after the first earthquake is Aj ~ 0.14 TECU at the equilibrium value of / ~ 5 
TECU. Assuming that on the order of magnitude Aj/I ~ u/C for a maximum gas 
velocity in the wave, we have an estimate of u exp ~ 30 m/s. For the sinusoidal form of 
the bipolar pulse with a duration T exp ~ 350 s, we find the experimental value of gas 
displacement along the direction of wave propagation: x+ exp = u exp T exp /tt ~ 3 km. 
Using the relation x+ exp ~ X+ it is possible to estimate the vertical displacement £o 
of the terrestrial surface in the epicentral zone of the earthquake. For this purpose, it 
seems reasonable to introduce the assumption about a short velocity impulse which 
models the main earthquake shock, rj <C 1, although the numerical value of x+ is 
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virtually independent of the value of 77 (|l4|). In view of the above considerations, we 
then obtain: 

Co ~ x+ exp ^RismP/{AL). (15) 

The uncertainty in the determination of £0 is caused both by the uncertainty of 
the true values of the quantities involved in this relation and by the limitations of 
the acoustic signal generation model itself. In this case, of the greatest importance 
is the dependence of £0 on the value of the acoustic parameter A, containing the 
atmospheric density p at the effective height z = z e ft which we have introduced 
artificially. The employment of the MSISE90 atmospheric model (Hedin, 1991), 
calculated for the location and time of the first earthquake, gives the values of 
£0 ~ 60,40 and 25 cm, with the values of z e jf = 350,400 and 450 km, respectively. 
In all cases it was assumed that R\ ~ &h0^f2 : K km, (3 = 20°, L = 30 km. The 
possibility of a vertical displacement of the terrestrial surface in the epicentral zone 
by several tens of centimeters seems real. In particular, (Calais and Minster, 1995) 
give the value of £0 ~ 40 cm at the epicenter of the Mw=6.7 Northridge earthquake 
(California, 1994). 

In view of the demonstration character of the above calculation, we have inten- 
tionally excluded from consideration the effects associated with the inclination of the 
magnetic field lines, and with the possible presence of strong winds at ionospheric 
heights. The presence of a magnetic field modifies the picture of transfer of move- 
ments from the neutral gas to the electron component of the ionosphere. Since the 
magnetic field is not entrained by the neutral gas, the field lines can be considered 
fixed. In this case the acceptable approximation would be the one, in which the 
electron component travels only along magnetic field lines with the velocity u-cosip, 
where is the angle between the magnetic field vector and the velocity vector of the 
neutral gas. Therefore, the quantity u ■ cosip must be involved in lieu of the quantity 
u in the expression Aj/I ~ u/C that was used above. 
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Let us estimate the value of ip for Turkey earthquakes. In the examples under 
discussion (Table 3), the horizontal projection K of the full wave vector K% is virtu- 
ally collinear to the horizontal component of the magnetic field. The wave vectors 
Kt, in turn, make with the vertical the angles from (3\ ~ 70° to @2 ~ 50°. 

Since the magnetic dip in the middle of Turkey is about 60° (the angle is mea- 
sured from the horizontal plane), tp = 20°-40°, and the value of costp ~ 0.94-0.77 
differs little from 1. It should be remembered, however, that taking into account 
this factor can be very important in the analysis of the complete picture of TEC 
disturbances above the earthquake or explosion source (see, for example, Calais et 
al., 1998). 

The presence of the zonal and meridional winds at ionospheric heights leads to a 
displacement and deformation of the wave front, and hence give rise to a dependence 
of the acoustic wave intensity on the propagation direction. The decisive role in 
this case is played by the wind velocity gradient. This factor can be taken into 
account within the framework of the ray theory. However, a corresponding model 
calculation would be worthwhile in the analysis of experimental data obtained for 
a set of subionospheric points surrounding the acoustic wave source. In the present 
situation, however, where the number of subionospheric points used in the analysis 
is too small, and the uncertainty of the parameters of the acoustic emitter itself 
is too large, the solution of such an unwieldy problem would be an overrun of the 
accuracy which is pursued by the above computational scheme. 



As follows from the expressions (13), maximum values of displacements and of the 
velocity of the neutral atmospheric species are attained directly above the earthquake 
epicenter. The signal duration is minimal, and does not seem to exceed a few tens of 
seconds at ionospheric heights. Since the wave vector of the disturbance is directed 
predominantly upward, the method of oblique-incidence ionospheric sounding in this 
case is a technique of choice for determining the waveform. 
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6 Conclusion 



In this paper we have investigated the form and dynamics of shock-acoustic waves 
generated during earthquakes. We have developed a method of determining the 
SAW parameters using GPS arrays whose elements can be chosen out of a large set 
of the global network GPS stations. Unlike existing radio techniques, the proposed 
method estimates the SAW parameters without a priori information about the lo- 
cation and time of the earthquake. The implementation of the method is illustrated 
by an analysis of ionospheric effects of the earthquakes in Turkey (August 17, and 
November 12, 1999), in Southern Sumatera (June 4, 2000), and off the coast of 
Central America (January 13, 2001). 

It was found that, in spite of the difference of the earthquake characteristics, the 
local time, the season, and the level of geomagnetic disturbance, for four earthquakes 
the time period of the ionospheric response is 180-390 s, and the amplitude exceeds 
by a factor of two as a minimum the standard deviation of background fluctuations in 
total electron content in this range of periods under quiet and moderate geomagnetic 
conditions. 

As has been pointed out in the Introduction, some investigators report markedly 
differing values of the SAW propagation velocity — by as much as several thousands 
m/s, which is beyond the values of the sound velocity at the SAW propagation 
heights in the atmosphere. The method proposed in this paper opens up a possibility 
of determining the angular characteristics of the wave vector K t and, accordingly, 
of estimating V t . According to our data (Table 3), the elevation of the SAW wave 
vector varied within 20-44° , and the phase velocity of the SAW varied from 1100 to 
1300 m/s. We determine the phase velocity of the equal TEC line at the height of the 
ionospheric F-region maximum which makes the main contribution to variations of 
the TEC between the receiver and the GPS satellite, and corresponds to the region 
of maximum sensitivity of the method. Since Vt approaches the sound velocity at 
these heights (Li et al., 1994) this makes it possible to identify the sound origin of the 
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TEC disturbance. The SAW source location, calculated without taking into account 
the refraction corrections, approximately corresponds to the earthquake epicenter. 
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TURKEY EARTHQUAKES 
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EARTHQUAKE OFF COAST OF CENTRAL AMERICA 
AND SOUTHERN SUMATERA 




Fig. 1: Experimental geometry during the earthquakes in Turkey - a). Crosses 
show the positions of the earthquake epicenters. Solid curves represent the tra- 
jectories of the subionospheric points for each GPS satellite at the height /ij=400 
km. Dark diamonds along the trajectories correspond to the coordinates of the 
subionospheric points at time t p of a maximum deviation of the TEC. Heavy dots 
and large lettering show the location and the names of the GPS stations, while 
lower-case letters along the trajectories refer to station names and PRN numbers 
of the GPS satellites for these trajectories. 

Asterisks mark the source location at km altitude inferred from the data from 
the GPS arrays. Numbers at the asterisks correspond to the respective day num- 
bers. Straight dashed lines that connect the expected source and the subiono- 
spheric point represent the horizontal projection of the wave vector K t . The 
scaling of the coordinate axes is chosen from considerations of an approximate 
equality of the linear dimensions along the latitude and longitude. 
The same, but for the earthquakes off the coast of Central America on January 
13, 2001 - b), and in Southern Sumatera on June 4, 2000- c). 
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Fig. 2: Time dependencies of slant TEC I(t) at one of the three sites of the 
GPS-array in the area of the earthquakes on August 17 and November 12, 1999, 
on the day of earthquake (heavy curve), and one day before and one day after 
the event (thin curves) — panels a, d. Panels b, e — for the same days but the 
TEC variations AJ(£) with the linear trend removed and with a smoothing with 
a time window of 5 min; panels c, f — variations of the frequency Doppler shift 
F(t) 'reduced' to the sounding signal frequency of 136 MHz, for three sites of the 
arrays, days of earthquake. All panels show day numbers, GPS station names, 
and PRN numbers of the GPS satellites. 



EARTHQUAKES OFF COAST OF CENTRAL AMERICA 
AND SOUTHERN SUMATERA 
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Fig. 3: Same as in Fig. 2, but for the earthquakes off the coast of Central America, 
January 13, 2001 - at the left, and in Southern Sumatera, June 4, 2000 - at the 
right. 
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Fig. 4: a. Model time dependence (from top to bottom) of the vertical displace- 
ment £, the velocity U = d^/dt and the acceleration a = dll/dt of the terrestrial 
surface in the epicentral zone of the earthquake. 

b. Acoustic signals u(t, (3) in the far-field radiation zone of the piston 2L = 60 km 
in diameter. The piston's velocity has the form of a rectangular impulse of a 
duration At = 10 s. The signals correspond to the expression (12) under the 
assumption of a homogeneous atmosphere; A = 1, and C = Co ~ 0.34 km/s. 
The zenith angles are (3 = 10°, 20°, 40°, 60°, and 90°. The abscissa axis indicates 
the time r in seconds, and the axis of ordinates indicates the gas velocity u. The 
amplitude of the strongest signal (3 = 10° is taken to be unity. 




Table 1: General information about earthquakes 



N 


Epicenter 


Data 


to 
(UT) 


Depth 
(km) 


Magnitude 
mb Ms Mw 


DST 
(nT) 


1 


40.70 N, 29.99 E 


Aug. 17 1999 (DAY 229) 


00 


01:39 


17 


6.3 7.8 7.4 


-14 


2 


40.79 N, 31.11 E 


Nov. 12 1999 (DAY 316) 


16 


57:20 


10 


6.5 7.5 7.1 


-44 


3 


4.72 S, 102.1 E 


Jun. 4 2000 (DAY 156) 


16 


28:26 


33 


6.8 8.0 7.7 


+8 


4 


12.83 N, 88.79 W 


Jan. 13 2001 (DAY 013) 


17 


33:29 


39 


- - 7.6 


+4 



Table 2: General information about earthquakes 



N 


Sites 


Geograph. 
latitude 


Geograph. 
longitude 


N 


Sites 


Geograph. 
latitude 


Geograph. 
longitude 


1 


ANKR 


39.887 


32.759 


9 


TEGU 


14.090 


-87.206 


2 


BSHM 


32.778 


35.022 


10 


SSIA 


13.697 


-89.117 


3 


GILB 


32.479 


35.416 


11 


MANA 


12.149 


-86.2489 


4 


KABR 


33.022 


35.145 


12 


GUAT 


14.590 


-90.520 


5 


KATZ 


32.995 


35.688 


13 


ESTI 


13.100 


-86.3621 


6 


NICO 


35.141 


33.396 


14 


SAMP 


3.622 


98.715 


7 


RAMO 


30.597 


34.763 


15 


NTUS 


1.346 


103.680 


8 


TELA 


32.067 


34.780 


16 


BAKO 


-6.491 


106.849 



Table 3: The parameters of shock-acoustic waves 



No. 


Sites 


tp 
(UT) 


At, 
sec. 


T, 
sec. 


At, 
TECU 


A F , 
Hz 


0,° 


a,° 


v h , 

m/s 


V u 
m/s 


V a , 
m/s 


o 

<Pw> 

° 










1999 Aug. 17; 


^0=00:01:39 UT 










1 


KABR 
BSHM 
KATZ 


00:21:04 
00:21:19 


1165 
1180 


360 
300 


0.15 
0.14 


0.037 
0.049 


24.9 


155 


1296 


1174 


873 
862 


39.6 
26.4 


2 


GILB 

KATZ 

BSHM 


00:21:54 


1215 


390 


0.12 


0.043 


26.1 


154 


1307 


1174 


868 


39.1 
25.9 


3 


KATZ 
KABR 
GILB 


00:21:23 


1184 


360 


0.19 


0.05 


25.1 


155 


1303 


1179 


854 


39.5 
26.4 


4 


TELA 
BSHM 
GILB 


00:22:14 


1235 


360 


0.1 


0.026 


19.9 


161 


1238 


1164 


894 


41.5 
26.1 


5 


E 






354 


0.14 


0.04 


24.1 


156 


1286 


1173 


870 


39.9 
26.2 


1999 Nov. 12; t = 16:57:20 UT 


6 


KABR 
BSHM 
GILB 


17:12:58 
17:13:21 


938 
961 


180 
180 


0.06 
0.07 


0.027 
0.021 


29.5 


194 


1285 


1119 


807 
787 


42.2 
30.1 


7 


TELA 
KABR 
GILB 


17:14:09 


1009 


210 


0.097 


0.021 


43.8 


180 


1487 


1073 


838 


39.4 
28.6 


8 


GILB 

BSHM 

TELA 


17:13:39 


979 


210 


0.09 


0.022 


39.6 


184 


1663 


1280 


817 


40.1 
29.1 


9 


E 






195 


0.079 


0.023 


37.6 


186 


1478 


1157 


812 


40.5 
29.2 










2000 Jun. 4; 


t =16:28:26 UT 










10 


SAMP 
NTUS 
BAKO 


16:46:12 
16:42:16 
16:49:16 


1092 

856 

1276 


270 
270 
240 


0.1 
0.5 
0.06 


0.04 
0.09 
0.02 










558 
503 
642 




2001 Jan. 13; i = 17:33:29 UT 


11 


TEGU 
ESTI 

MANA 


17:45:11 
17:46:30 
17:47:21 


851 
930 
981 


240 
210 
270 


0.09 
0.09 
0.2 


0.03 
0.03 
0.05 










591 
488 
440 





